Aminet 1/Aminet - June 1993 [Walnut Creek].iso/usenet/sources/volume90/aplictns/listplot/part03

Internet Message Format | 1990-08-23 | 28.1 KB

Path: abcfd20.larc.nasa.gov!amiga-request
From: [e-mail] (Amiga Sources/Binaries Moderator)
Subject: v90i244: ListPlot - 2d plotting program, Part03/03
Reply-To: Anthony M. Richardson <[e-mail]>
Newsgroups: comp.sources.amiga
Message-ID: <comp.sources.amiga:[e-mail]>
Date: 23 Aug 90 01:03:49 GMT
Approved: [e-mail] (Tad Guy)
X-Mail-Submissions-To: [e-mail]
X-Post-Discussions-To: comp.sys.amiga

Submitted-by: Anthony M. Richardson <[e-mail]>
Posting-number: Volume 90, Issue 244
Archive-name: applications/listplot/part03

#!/bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 3 (of 3)."
# Contents:  Csrc/plotting.c
# Wrapped by tadguy@abcfd20 on Wed Aug 22 21:03:22 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Csrc/plotting.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Csrc/plotting.c'\"
else
echo shar: Extracting \"'Csrc/plotting.c'\" \(25452 characters\)
sed "s/^X//" >'Csrc/plotting.c' <<'END_OF_FILE'
X/*
X * plotting.c
X *
X * plotting routines.
X *
X */
X#include <stdio.h>
X#include <errno.h>
Xextern int errno;
X#include <assert.h>
X#include <string.h>
X#include <ctype.h>
X#include <signal.h>
X#include "datatypes.h"
X
Xvoid    ErrorExit();
X
X#ifndef    min
X#    define    min(A,B)    A<B? A : B
X#endif
X#ifndef    max
X#    define    max(A,B)    A>B? A : B
X#endif
X
Xextern bool    GraphicsInProgress;
X
X/* GoldenRatio = (1+Sqrt(5))/2    */
X#define    GOLDENRATIO    1.6180339887499
X#define ASPECTRATIO    1.0/GOLDENRATIO
Xextern VALUE    V_AngularUnit;
Xextern VALUE    V_AnnotationScale;
Xextern VALUE    V_AspectRatio;
Xextern VALUE    V_Boxed;
Xextern VALUE    V_Domain;
Xextern VALUE    V_Gridding;
Xextern VALUE    V_LabelScale;
Xextern VALUE    V_LineColor;
Xextern VALUE    V_LineStyle;
Xextern VALUE    V_Orientation;
Xextern VALUE    V_Origin;
Xextern VALUE    V_PlotColor;
Xextern VALUE    V_PlotDevice;
Xextern VALUE    V_PlotJoined;
Xextern VALUE    V_PlotPoints;
Xextern VALUE    V_PlotTitle;
Xextern VALUE    V_PlotType;
Xextern VALUE    V_PointScale;
Xextern VALUE    V_PointSymbol;
Xextern VALUE    V_PolarVariable;
Xextern VALUE    V_Range;
Xextern VALUE    V_SubPages;
Xextern VALUE    V_SupplyAbscissa;
Xextern VALUE    V_TitleScale;
Xextern VALUE    V_UseInputFile;
Xextern VALUE    V_Verbose;
Xextern VALUE    V_ViewPort;
Xextern VALUE    V_XLabel;
Xextern VALUE    V_XTick;
Xextern VALUE    V_YLabel;
Xextern VALUE    V_YTick;
X
Xchar    *DeviceTypes[] = {
X    "amiga",
X    "printer",
X    "iff",
X    "hp",
X    "aegis",
X    "postscript",
X    (char *)NULL
X};
X
Xtypedef struct line_style {
X        int    ls_type;
X        int    ls_nelms;
X        int    *ls_space;
X        int    *ls_mark;
X    } LINESTYLE;
XLIST    LineStyleList;
XLINESTYLE    *LineStyles; int NLineStyles;
Xint        *PointCodes, NPointCodes;
Xint    space0 = 0, mark0 = 0, space1 = 1500, mark1 = 1500;
X
XRECT    UserRect;
XRECT    ViewPort = { 0.1,0.1, 0.9,0.9 };
XRECT    WorldRect = {0.0,0.0, 1.0,1.0};
X
X
Xvoid
XListPlot(Data, NTuples, TupleSize)
XFLOAT    **Data;
Xint    NTuples, TupleSize;
X{
X    void InitializeDevice();
X    void InitializePlottingEnv();
X#ifdef    ANSI_C
X    void Plot2D(FLOAT **, int , int );
X    void PolarPlot(FLOAT **, int , int );
X#else
X    void Plot2D();
X    void PolarPlot();
X#endif
X
X    InitializeDevice();
X    GraphicsInProgress = TRUE;
X
X    InitializePlottingEnv(Data, NTuples, TupleSize);
X
X    if (Vstrcmp(V_PlotType,"polar")) {
X        PolarPlot(Data, NTuples, TupleSize);
X    } else {
X        Plot2D(Data, NTuples, TupleSize);
X    }
X    plend();
X}
X
X
Xvoid
XInitializeDevice()
X{
X    register int    i;
X    register char    *DeviceStr;
X
X
X    plorient(Vstrcmp(V_Orientation, "portrait")? 1 : 0);
X    plselect(stdout);    /* write to stdout */
X
X
X    DeviceStr = V_PlotDevice.val_u.udt_string;
X    for (i=0; DeviceTypes[i]; i++)
X        if (strcmp(DeviceStr, DeviceTypes[i]) == 0) {
X            plbeg(
X                i+1,    /* device code    */
X                (int)V_SubPages.val_u.udt_interval.int_lo, /* nx */
X                (int)V_SubPages.val_u.udt_interval.int_hi /* ny */
X            );
X            return;
X        }
X    fprintf(stderr,"(InitializeDevice) Unknown device type: %s\n", DeviceStr);
X    ErrorExit();
X}
X
Xvoid
XInitializePlottingEnv(Data, NTuples, TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X{
X    int    i, j;
X    FLOAT    *Ptr;
X    void    InitLineStyles();
X    void    InitPointCodes();
X
X    /* viewport    */
X    ViewPort = VtoRect(V_ViewPort);
X
X    /* user window        */
X    UserRect.XMin = MAX_FLOAT;
X    UserRect.YMin = MAX_FLOAT;
X    UserRect.XMax = -MAX_FLOAT;
X    UserRect.YMax = -MAX_FLOAT;
X
X    if (VtoBoolean(V_SupplyAbscissa)) {
X        UserRect.XMin = 0.0;
X        UserRect.XMax = (FLOAT)(NTuples-1);
X    } else {
X        for (i=0; i<NTuples; i++) {
X            if (Data[i][0] < UserRect.XMin)
X                UserRect.XMin = Data[i][0];
X            if (Data[i][0] > UserRect.XMax)
X                UserRect.XMax = Data[i][0];
X        }
X    }
X    for (i=0; i<NTuples; i++)
X        for (j=TupleSize-1,Ptr= &Data[i][1]; j>0; --j,Ptr++) {
X            if (*Ptr < UserRect.YMin)
X                UserRect.YMin = *Ptr;
X            if (*Ptr > UserRect.YMax)
X                UserRect.YMax = *Ptr;
X        }
X
X
X    /* font scaling        */
X
X    /* plot type        */
X
X    /* grid            */
X
X    /* line styles        */
X    InitLineStyles();
X
X
X    /* line colors        */
X
X    /* point symbols    */
X    InitPointCodes();
X
X    /* window type        */
X
X    /* X label        */
X
X    /* Y label        */
X
X    /* Plot title        */
X
X    /* aspect ratio        */
X}
X
X#define    isMark(C)    (((C) == 'M' || (C) == 'm')? TRUE : FALSE)
X#define    isSpace(C)    (((C) == 'S' || (C) == 's')? TRUE : FALSE)
X#define    MarkLength(C)    (C) == 'M'? 1000 : 250
X#define    SpaceLength(C)    (C) == 'S'? 1000 : 250
X#define    MAX_MARKS    32
X
Xvoid
XDecodeLineStyle(StyleStr, Mark, Space, N)
Xchar    *StyleStr;
Xint    **Mark;
Xint    **Space;
Xint    *N;
X{
X    register char    *S;
X    int    macc, sacc;
X    int    i, j, k;
X    int    mark[MAX_MARKS];
X    int    space[MAX_MARKS];
X    bool    ProcessingMark;
X
X    /*
X     * 'S' = 1000   micron space
X     * 's' =  250   micron space
X     * 'M' = 1000     micron line
X     * 'm' =  250    micron line
X     */
X    ProcessingMark = TRUE; /* to satisfy some compilers */
X    if (isMark(*StyleStr))
X        ProcessingMark = TRUE;
X    else if (isSpace(*StyleStr))
X        ProcessingMark = FALSE;
X    else {
X        fprintf(stderr,"(DecodeLineStyle) Invalid Line style: \"%s\"\n", StyleStr);
X        ErrorExit();
X    }
X    memset((char *)mark, 0, MAX_MARKS * sizeof(int));
X    memset((char *)space, 0, MAX_MARKS * sizeof(int));
X    for (macc=0,sacc=0,i=0,j=0,S=StyleStr ; *S && i<MAX_MARKS && j<MAX_MARKS; ++S) {
X        if (isMark(*S) && ProcessingMark) {
X            macc += MarkLength(*S);
X        } else if (isMark(*S) && !ProcessingMark) {
X            space[j++] = sacc;
X            macc = MarkLength(*S);
X            ProcessingMark = TRUE;
X        } else if (isSpace(*S) && ProcessingMark) {
X            mark[i++] = macc;
X            sacc = SpaceLength(*S);
X            ProcessingMark = FALSE;
X        } else if (isSpace(*S) && !ProcessingMark) {
X            sacc += SpaceLength(*S);
X        } else if (isspace(*S)) {
X            continue;
X        } else {
X            fprintf(stderr,"(DecodeLineStyle) Invalid line style code '%c' in \"%s\"\n", *S, StyleStr);
X            ErrorExit();
X        }
X    }
X    if (ProcessingMark)
X        mark[i++] = macc;
X    else
X        space[j++] = sacc;
X
X    k = max(i,j);
X
X    if (((*Mark) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
X        perror("(DecodeLineStyle) ");
X        ErrorExit();
X    }
X    if (((*Space) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
X        perror("(DecodeLineStyle) ");
X        ErrorExit();
X    }
X    for (i=0; i<k; i++) {
X        (*Mark)[i] = mark[i];
X        (*Space)[i] = space[i];
X    }
X    *N = k;
X}
X
Xvoid
XInitLineStyles()
X{
X    int        i;
X    LATOM        *A;
X    LINESTYLE    *LS;
X    int        *Mark, *Space;
X    int        N;
X#ifdef    ANSI_C
X    LINESTYLE    *LSAlloc(int , int *, int *);
X#else
X    LINESTYLE    *LSAlloc();
X#endif
X
X    LineStyleList.list_n = 0;
X    LineStyleList.list_head = (LATOM *)NULL;
X    LineStyleList.list_tail = (LATOM *)NULL;
X
X    LS = LSAlloc(0, (int *)NULL, (int *)NULL);
X    Append(&LineStyleList, (generic *)LS);
X
X    for (i=1,A=V_LineStyle.val_u.udt_set.list_head; A; A=A->la_next,i++) {
X        DecodeLineStyle(
X            (char *)A->la_p,
X            &Mark,
X            &Space,
X            &N
X        );
X        LS = LSAlloc(N, Mark, Space);
X        Append(&LineStyleList, (generic *)LS);
X    }
X
X    NLineStyles = i;
X    if ((LineStyles = (LINESTYLE *)calloc(NLineStyles, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
X        perror("(InitLineStyles) ");
X        ErrorExit();
X    }
X    for (i=0,A=LineStyleList.list_head; A; A=A->la_next,i++) {
X        LineStyles[i].ls_mark = ((LINESTYLE *)(A->la_p))->ls_mark;
X        LineStyles[i].ls_space = ((LINESTYLE *)(A->la_p))->ls_space;
X        LineStyles[i].ls_nelms = ((LINESTYLE *)(A->la_p))->ls_nelms;
X    }
X}
X
X
XLINESTYLE    *
XLSAlloc(n, Mark, Space)
Xint    n;
Xint    *Mark, *Space;
X{
X    register LINESTYLE    *LS;
X
X    if ((LS = (LINESTYLE *)calloc(1, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
X        perror("(LSAlloc) ");
X        assert(LS != (LINESTYLE *)NULL);
X        exit(errno);
X    }
X    LS->ls_mark = Mark;
X    LS->ls_space = Space;
X    LS->ls_nelms = n;
X    return(LS);
X}
X
Xvoid
XInitPointCodes()
X{
X    register LATOM    *A;
X    register int    i;
X    extern int    *PointCodes, NPointCodes;
X    extern VALUE    V_PointSymbol;
X
X    if (!isSet(V_PointSymbol)) {
X        NPointCodes = 0;
X        return;
X    }
X
X    /* count # codes    */
X    for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++)
X        ;
X
X    NPointCodes = i;
X    if (NPointCodes == 0)
X        return;
X
X    if ((PointCodes = (int *)calloc(NPointCodes, sizeof(int))) == (int *)NULL) {
X        perror("(InitPointCodes) ");
X        assert(PointCodes != (int *)NULL);
X        ErrorExit();
X    }
X    for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++) {
X        PointCodes[i] = atoi((char *)A->la_p);
X    }
X}
X
X
Xvoid
XPlot2D(Data, NTuples, TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X{
X    int    i, j;
X    FLOAT    *X, *Y;
X    FLOAT    x, y;
X    FLOAT    x0, y0;
X    FLOAT    xtick, ytick;
X    FLOAT    Xmin,Xmax, Ymin,Ymax;
X    FLOAT    MedianX, StdDevX;
X    FLOAT    MedianY, StdDevY;
X    int    nxsub, nysub;
X    bool    AutoRange, AutoDomain;
X    bool    LogX, LogY;
X    char    Xopt[16], Yopt[16];
X#ifdef    ANSI_C
X    void    StatisticsOfX(
X            FLOAT **,
X            int,
X            int ,
X            FLOAT *,
X            FLOAT *
X        );
X    void    StatisticsOfY(
X            FLOAT **,
X            int,
X            int ,
X            FLOAT *,
X            FLOAT *
X        );
X    double    log10(double);
X#else
X    void    StatisticsOfX();
X    void    StatisticsOfY();
X    double    log10();
X#endif
X
X    pladv(0);
X    ViewPort = VtoRect(V_ViewPort);
X
X    if (isString(V_AspectRatio)) {
X        /* automatic --> STRETCH */
X        plvpor(ViewPort.XMin, ViewPort.XMax, ViewPort.YMin, ViewPort.YMax);
X/*debug
X        plvsta();
Xdebug*/
X    } else {
X        assert(isDbl(V_AspectRatio));
X        plvasp(VtoDbl(V_AspectRatio));
X    }
X
X    LogX = (Vstrcmp(V_PlotType, "loglin") || Vstrcmp(V_PlotType,"loglog"))?
X            TRUE : FALSE;
X    LogY = (Vstrcmp(V_PlotType, "linlog") || Vstrcmp(V_PlotType,"loglog"))?
X            TRUE : FALSE;
X
X    /* window bounds */
X    AutoDomain = FALSE;
X    if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
X        /* all points    */
X        Xmin = UserRect.XMin;
X        Xmax = UserRect.XMax;
X    } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
X        /* automatic     */
X        AutoDomain = TRUE;
X        StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
X        x = MedianX - 3.0*StdDevX;
X        Xmin = UserRect.XMin < x? x : UserRect.XMin;
X        x = MedianX + 3.0*StdDevX;
X        Xmax = UserRect.XMax > x? x : UserRect.XMax;
X    } else {
X        assert(isInterval(V_Domain));
X        Xmin = V_Domain.val_u.udt_interval.int_lo;
X        Xmax = V_Domain.val_u.udt_interval.int_hi;
X    }
X    if (LogX) {
X        if (Xmin == 0.0 || Xmax == 0.0) {
X            /* error */
X            fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
X            ErrorExit();
X        }
X        Xmin = log10((double)Xmin);
X        Xmax = log10((double)Xmax);
X    }
X
X    AutoRange = FALSE;
X    if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
X        /* all */
X        Ymin = UserRect.YMin;
X        Ymax = UserRect.YMax;
X    } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
X        /* automatic     */
X        AutoRange = TRUE;
X        StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
X        y = MedianY - 3.0*StdDevY;
X        Ymin = UserRect.YMin < y? y : UserRect.YMin;
X        y = MedianY + 3.0*StdDevY;
X        Ymax = UserRect.YMax > y? y : UserRect.YMax;
X    } else {
X        assert(isInterval(V_Range));
X        Ymin = V_Range.val_u.udt_interval.int_lo;
X        Ymax = V_Range.val_u.udt_interval.int_hi;
X    }
X    if (LogY) {
X        if (Ymin == 0.0 || Ymax == 0.0) {
X            /* error */
X            fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
X            ErrorExit();
X        }
X        Ymin = log10((double)Ymin);
X        Ymax = log10((double)Ymax);
X    }
X    plwind(
X        Xmin,
X        Xmax,
X        Ymin,
X        Ymax
X    );
X    plschr(0., VtoDbl(V_AnnotationScale));
X
X
X    /* X and Y ticks */
X    if (isString(V_XTick)) {
X        /* automatic */
X        xtick = 0.;
X        nxsub = 0;
X    } else {
X        xtick = V_XTick.val_u.udt_interval.int_lo;
X        if (LogX) xtick = log10((double)xtick);
X        nxsub = V_XTick.val_u.udt_interval.int_hi;
X    }
X    if (isString(V_YTick)) {
X        /* automatic */
X        ytick = 0.;
X        nysub = 0;
X    } else {
X        ytick = V_YTick.val_u.udt_interval.int_lo;
X        if (LogY) ytick = log10((double)ytick);
X        nysub = V_YTick.val_u.udt_interval.int_hi;
X    }
X
X    /* Axis options    */
X    Xopt[0] = (char)NULL;
X    Yopt[0] = (char)NULL;
X    if (isBoolean(V_Boxed) && VtoBoolean(V_Boxed)) {
X        strcat(Xopt, "bc");
X        strcat(Yopt, "bc");
X    } else if (isString(V_Boxed)) {
X        if (Vstrchr(V_Boxed, 'b'))
X            strcat(Xopt, "b");
X        if (Vstrchr(V_Boxed, 't'))
X            strcat(Xopt, "c");
X        if (Vstrchr(V_Boxed, 'r'))
X            strcat(Yopt, "c");
X        if (Vstrchr(V_Boxed, 'l'))
X            strcat(Yopt, "b");
X    } else {
X        ;
X    }
X    if (LogX)
X        strcat(Xopt, "l");
X    if (LogY)
X        strcat(Yopt, "l");
X    strcat(Xopt, "nst");
X    strcat(Yopt, "nstv");
X
X    if (isString(V_Origin) && Vstrcmp(V_Origin, "Median")) {
X        if (AutoDomain == FALSE) {
X            StatisticsOfX(Data,NTuples,TupleSize,&MedianX,&StdDevX);
X        }
X        if (AutoRange == FALSE) {
X            StatisticsOfY(Data,NTuples,TupleSize,&MedianY,&StdDevY);
X        }
X        plaxes(
X            MedianX,MedianY,
X            Xopt, xtick, nxsub,
X            Yopt, ytick,nysub
X        );
X    } else if (isInterval(V_Origin)) {
X        x0 = V_Origin.val_u.udt_interval.int_lo;
X        y0 = V_Origin.val_u.udt_interval.int_hi;
X        plaxes(
X            x0,y0,
X            Xopt, xtick, nxsub,
X            Yopt, ytick,nysub
X        );
X    } else if (isString(V_Origin) && Vstrcmp(V_Origin, "Automatic")) {
X        plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
X    } else {
X        plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
X    }
X
X    /* gridding    */
X    if (VtoBoolean(V_Gridding)) {
X        /* gridding on */
X        plstyl(1, &mark1, &space1);
X        plbox("g", xtick, nxsub, "g", ytick,nysub);
X        plstyl(0, &mark0, &space0);
X    }
X
X    /* plot curves    */
X    if ((X = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X        perror("(Plot2D) ");
X        ErrorExit();
X    }
X    if ((Y = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X        perror("(Plot2D) ");
X        ErrorExit();
X    }
X    if (VtoBoolean(V_SupplyAbscissa)) {
X        for (i=0; i<NTuples; i++)
X            X[i] = LogX? log10((double)(i+1)) : (FLOAT)i;
X    } else {
X        for (i=0; i<NTuples; i++)
X            X[i] = LogX? log10((double)Data[i][0]) : Data[i][0];
X    }
X    for (i=1; i<TupleSize; i++) {
X        if (NLineStyles > 0) {
X            /* use user-supplied patterns */
X            plstyl(
X                LineStyles[(i-1)%NLineStyles].ls_nelms,
X                LineStyles[(i-1)%NLineStyles].ls_mark,
X                LineStyles[(i-1)%NLineStyles].ls_space
X            );
X        } else {
X            pllsty((i-1)%8 + 1);
X        }
X        plcol(i);
X        for (j=0; j<NTuples; j++)
X            Y[j] = LogY? log10((double)Data[j][i]) : Data[j][i];
X
X        if (VtoBoolean(V_PlotJoined))
X            plline(NTuples, X, Y);
X
X        if (VtoBoolean(V_PlotPoints)) {
X            plschr(0., VtoDbl(V_PointScale));
X            pllsty(1);
X            if (NPointCodes > 0) {
X                plpoin(NTuples, X, Y, PointCodes[(i-1)%NPointCodes]);
X            } else {
X                plpoin(NTuples, X, Y, i);
X            }
X        }
X    }
X
X
X    /* restore line styles etc... */
X    pllsty(1);
X    plcol(1);
X    plschr(0., VtoDbl(V_LabelScale));
X    pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
X    plschr(0., VtoDbl(V_TitleScale));
X    pllab("","", VtoString(V_PlotTitle));
X}
X
X#define    DTR(D)    (double)((D)*3.1415927/180.0)
X
Xvoid
XPolarPlot(Data, NTuples, TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X{
X    int    i, j;
X    FLOAT    *A, *R;
X    FLOAT    x, y;
X    FLOAT    xtick, ytick;
X    FLOAT    atick, rtick;
X    FLOAT    Xmin,Xmax, Ymin,Ymax;
X    FLOAT    r, rmax, ang, amax;
X    FLOAT    MedianX, StdDevX;
X    FLOAT    MedianY, StdDevY;
X    FLOAT    da;
X    FLOAT    x0,y0, xn,yn;
X    double    pi;
X    char    s[32];
X    int    nxsub, nysub, nrsub, nasub;
X    bool    InRadians, XisAngular;
X    bool    AutoDomain, AutoRange;
X#ifdef    ANSI_C
X    void    StatisticsOfX(
X            FLOAT **,
X            int ,
X            int ,
X            FLOAT *,
X            FLOAT *
X        );
X    void    StatisticsOfY(
X            FLOAT **,
X            int ,
X            int ,
X            FLOAT *,
X            FLOAT *
X        );
X    double    log10(double); double fabs();
X    double    sqrt(), sin(), cos(), atan();
X#else
X    void    StatisticsOfX();
X    void    StatisticsOfY();
X    double    log10(); double fabs();
X    double    sqrt(), sin(), cos(), atan();
X#endif
X
X    pi = 4.0 * atan((double)1.0);
X    InRadians = Vstrcmp(V_AngularUnit, "degrees")? FALSE : TRUE;
X    XisAngular = Vstrcmp(V_PolarVariable,"angle")? TRUE : FALSE;
X
X    pladv(0);
X
X    if (isString(V_AspectRatio)) {
X        /* automatic --> STRETCH */
X        plvasp(1.0);
X    } else {
X        assert(isDbl(V_AspectRatio));
X        plvasp(VtoDbl(V_AspectRatio));
X    }
X
X    /* window bounds */
X    AutoDomain = FALSE;
X    if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
X        /* all points    */
X        Xmin = UserRect.XMin;
X        Xmax = UserRect.XMax;
X    } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
X        /* automatic     */
X        AutoDomain = TRUE;
X        StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
X        x = MedianX - 3.0*StdDevX;
X        Xmin = UserRect.XMin < x? x : UserRect.XMin;
X        x = MedianX + 3.0*StdDevX;
X        Xmax = UserRect.XMax > x? x : UserRect.XMax;
X    } else {
X        assert(isInterval(V_Domain));
X        Xmin = V_Domain.val_u.udt_interval.int_lo;
X        Xmax = V_Domain.val_u.udt_interval.int_hi;
X    }
X
X    AutoRange = FALSE;
X    if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
X        /* all */
X        Ymin = UserRect.YMin;
X        Ymax = UserRect.YMax;
X    } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
X        /* automatic     */
X        AutoRange = TRUE;
X        StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
X        y = MedianY - 3.0*StdDevY;
X        Ymin = UserRect.YMin < y? y : UserRect.YMin;
X        y = MedianY + 3.0*StdDevY;
X        Ymax = UserRect.YMax > y? y : UserRect.YMax;
X    } else {
X        assert(isInterval(V_Range));
X        Ymin = V_Range.val_u.udt_interval.int_lo;
X        Ymax = V_Range.val_u.udt_interval.int_hi;
X    }
X    rmax = XisAngular? Ymax : Xmax;
X    plwind(
X        -1.3*rmax,
X        1.3*rmax,
X        -1.3*rmax,
X        1.3*rmax
X    );
X
X
X    /* X and Y ticks */
X    if (isString(V_XTick)) {
X        /* automatic */
X        xtick = XisAngular? (InRadians? pi/6.0 : 30.0) : 0.25 * rmax;
X        nxsub = 4;
X    } else {
X        xtick = V_XTick.val_u.udt_interval.int_lo;
X        nxsub = V_XTick.val_u.udt_interval.int_hi;
X    }
X    if (isString(V_YTick)) {
X        /* automatic */
X        ytick = XisAngular? 0.25 * rmax : (InRadians? pi/6.0 : 30.0);
X        nysub = 4;
X    } else {
X        ytick = V_YTick.val_u.udt_interval.int_lo;
X        nysub = V_YTick.val_u.udt_interval.int_hi;
X    }
X    rtick = XisAngular? ytick : xtick;
X    atick = XisAngular? xtick : ytick;
X    nasub = XisAngular? nxsub : nysub;
X    nrsub = XisAngular? nysub : nxsub;
X    amax = InRadians? 2.0 * pi : 360.0;
X
X    /* Axis options    */
X    if (VtoBoolean(V_Gridding)) {
X        /* draw circles    */
X        plstyl(0, &mark1, space1);
X        /*  circular ticks    */
X        for (i=0,r=(rtick/nrsub); r<rmax; r += (rtick/nrsub),i++) {
X            da = i%nrsub? 0.05 * atick : 0.1 * atick;
X            for (ang=0; ang<amax; ang += atick) {
X                x0 = r * cos(InRadians? ang-da : DTR(ang-da));
X                y0 = r * sin(InRadians? ang-da : DTR(ang-da));
X                xn = r * cos(InRadians? ang+da : DTR(ang+da));
X                yn = r * sin(InRadians? ang+da : DTR(ang+da));
X                pljoin(x0,y0,xn,yn);
X            }
X        }
X#ifdef    CIRCULAR_TICKS
X        /* major circles    */
X        for (r=rtick; r<rmax; r += rtick) {
X            x0 = r; y0 = 0.0;
X            for (ang=1.0; ang<=360.0; ang++) {
X                xn = r * cos((double)(DTR(ang)));
X                yn = r * sin((double)(DTR(ang)));
X                pljoin(x0,y0,xn,yn);
X                x0 = xn; y0 = yn;
X            }
X        }
X#endif
X        /* draw spokes    */
X        for (i=0,ang=0.0; ang<amax; ang += (atick/nasub),i++) {
X            xn = rmax * cos(InRadians? ang : DTR(ang));
X            yn = rmax * sin(InRadians? ang : DTR(ang));
X            if (i%nasub == 0) {
X                pljoin((FLOAT)0.0,(FLOAT)0.0,xn,yn);
X            } else {
X                pljoin(xn,yn,1.05*xn,1.05*yn);
X            }
X        }
X    }
X    /* write angle labels    */
X    plschr(0., VtoDbl(V_AnnotationScale));
X    j = amax / atick;
X    for (ang=0.0,i=0; i<j; ang += atick,i++) {
X        if (InRadians) {
X            xn = rmax * cos((double)ang);
X            yn = rmax * sin((double)ang);
X            if (fabs(ang-(2.0*pi)) > 1.0e-2)
X                sprintf(s, "%.2f #gp", (float)ang/pi);
X        } else {
X            xn = rmax * cos((double)DTR(ang));
X            yn = rmax * sin((double)DTR(ang));
X            sprintf(s, "%d", (int)(ang+0.5));
X        }
X        if (xn >= 0)
X            plptex(xn,yn,xn,yn,-0.15, s);
X        else
X            plptex(xn,yn,-xn,-yn,1.15, s);
X    }
X    plstyl(0, &mark0, &space0);
X
X    /* plot curves    */
X    if ((A = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X        perror("(Plot2D) ");
X        ErrorExit();
X    }
X    if ((R = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X        perror("(Plot2D) ");
X        ErrorExit();
X    }
X    if (VtoBoolean(V_SupplyAbscissa)) {
X        for (i=0; i<NTuples; i++)
X            if (XisAngular)
X                A[i] = i*2.0*pi/(NTuples-1);
X            else
X                R[i] = i*rmax/(NTuples-1);
X    } else {
X        for (i=0; i<NTuples; i++)
X            if (XisAngular)
X                A[i] = InRadians? Data[i][0] :  DTR(Data[i][0]);
X            else
X                R[i] = Data[i][0];
X    }
X    for (i=1; i<TupleSize; i++) {
X        if (NLineStyles > 0) {
X            /* use user-supplied patterns */
X            plstyl(
X                LineStyles[(i-1)%NLineStyles].ls_nelms,
X                LineStyles[(i-1)%NLineStyles].ls_mark,
X                LineStyles[(i-1)%NLineStyles].ls_space
X            );
X        } else {
X            pllsty((i-1)%8);
X        }
X        plcol(i);
X
X        /* get and convert data if necessary    */
X        if (InRadians && !XisAngular)
X            for (j=0; j<NTuples; j++)
X                A[j] = Data[j][i];
X        else if ((!InRadians) && !XisAngular)
X            for (j=0; j<NTuples; j++)
X                A[j] = DTR(Data[j][i]);
X        else
X            for (j=0; j<NTuples; j++)
X                R[j] = Data[j][i];
X
X        /* convert to cartesian    A<->X  R<->Y*/
X        for (j=0; j<NTuples; j++) {
X            x0 = R[j] * cos(A[j]);
X            y0 = R[j] * sin(A[j]);
X            A[j] = x0;
X            R[j] = y0;
X        }
X        if (VtoBoolean(V_PlotJoined)) {
X            x0 = A[0];
X            y0 = R[0];
X            for (j=1; j<NTuples; j++) {
X                pljoin(A[j-1], R[j-1], A[j],R[j]);
X            }
X        }
X
X        if (VtoBoolean(V_PlotPoints)) {
X            plschr(0., VtoDbl(V_PointScale));
X            pllsty(1);
X            if (NPointCodes > 0) {
X                plpoin(NTuples, A, R, PointCodes[(i-1)%NPointCodes]);
X            } else {
X                plpoin(NTuples, A, R, i);
X            }
X        }
X    }
X
X    /* restore line styles etc... */
X    pllsty(1);
X    plcol(1);
X    plschr(0., VtoDbl(V_LabelScale));
X    pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
X    plschr(0., VtoDbl(V_TitleScale));
X    pllab("","", VtoString(V_PlotTitle));
X}
X
X
Xvoid
XStatisticsOfX(Data, NTuples,TupleSize, Median, StdDev)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
XFLOAT    *Median;
XFLOAT    *StdDev;
X{
X#ifdef    ANSI_C
X    FLOAT    MedianOfX();
X    FLOAT    StdDevOfX();
X/*debug
X    FLOAT    MedianOfX(FLOAT **Data,int NTuples,int TupleSize);
X    FLOAT    StdDevOfX(FLOAT **Data,int NTuples,int TupleSize);
Xdebug*/
X#else
X    FLOAT    MedianOfX();
X    FLOAT    StdDevOfX();
X#endif
X    *StdDev = StdDevOfX(Data,NTuples,TupleSize);
X    *Median = MedianOfX(Data,NTuples,TupleSize);
X}
X
X
Xvoid
XStatisticsOfY(Data, NTuples,TupleSize, Median, StdDev)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
XFLOAT    *Median;
XFLOAT    *StdDev;
X{
X#ifdef    ANSI_C
X    FLOAT    MedianOfY();
X    FLOAT    StdDevOfY();
X/*debug
X    FLOAT    MedianOfY(FLOAT **Data,int NTuples,int TupleSize);
X    FLOAT    StdDevOfY(FLOAT **Data,int NTuples,int TupleSize);
Xdebug*/
X#else
X    FLOAT    MedianOfY();
X    FLOAT    StdDevOfY();
X#endif
X    *StdDev = StdDevOfY(Data,NTuples,TupleSize);
X    *Median = MedianOfY(Data,NTuples,TupleSize);
X}
X
X
XFLOAT
XStdDevOfX(Data,NTuples,TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X{
X    int    i;
X    FLOAT    StdDev, x;
X    FLOAT    var, sum, mean;
X    double sqrt(), fabs();
X
X    if (VtoBoolean(V_SupplyAbscissa)) {
X        return((FLOAT)(NTuples*NTuples)/12.0);
X    }
X
X    sum = 0.0;
X    for (i=0; i<NTuples; i++)
X        sum += Data[i][0];
X    mean = sum / NTuples;
X
X    var = 0.0;
X    for (i=0; i<NTuples; i++) {
X        x = fabs((double)(Data[i][0] -  mean));
X        var += x * x;
X    }
X    var /= (NTuples-1);
X    StdDev = sqrt((double)var);
X    return(StdDev);
X}
X
X
XFLOAT
XStdDevOfY(Data,NTuples,TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X{
X    int    i,j,D;
X    FLOAT    StdDev, x;
X    FLOAT    var, sum, mean;
X    FLOAT    n;
X    double sqrt(), fabs();
X
X    D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
X    n = NTuples * (TupleSize-D);
X    sum = 0.0;
X    for (i=0; i<NTuples; i++)
X        for (j=D; j<TupleSize; j++) {
X            sum += Data[i][j];
X        }
X    mean = sum / n;
X
X    var = 0.0;
X    for (i=0; i<NTuples; i++)
X        for (j=D; j<TupleSize; j++) {
X            x = fabs((double)(Data[i][j] -  mean));
X            var += x * x;
X        }
X    var /= (n-1);
X    StdDev = sqrt((double)var);
X    return(StdDev);
X}
X
X#define    BIG    1.0e30
X#define    AFAC    1.5
X#define    AMP    1.5
X
XFLOAT
XMedianOfX(Data,NTuples,TupleSize)
XFLOAT    **Data;
Xint    NTuples;
Xint    TupleSize;
X/* an iterative computation of the median...
X * Code taken from Numerical Recipes..
X */
X{
X    int    np, nm, i;
X    FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
X    FLOAT    Median;
X    double    fabs();
X
X    if (VtoBoolean(V_SupplyAbscissa))
X        return((FLOAT)(0.5*NTuples));
X
X    /* first estimate    */
X    a = 0.5 * (Data[0][0] + Data[NTuples-1][0]);
X    eps = fabs((double)(Data[NTuples-1][0] - Data[0][0]));
X
X    am = -(ap=BIG);
X    for (;;) {
X        sum = sumx = 0.0;
X        np = nm = 0;        /* # point above and below median */
X        xm = -(xp = BIG);
X
X        for (i=0; i<NTuples; i++) {
X            xx = Data[i][0];
X            if (xx != a) {
X                if (xx > a) {
X                    ++np;
X                    if (xx < xp) xp = xx;
X                } else if (xx < a) {
X                    ++nm;
X                    if (xx > xm) xm = xx;
X                }
X                sum += dum=1.0/(eps+fabs((double)(xx-a)));
X                sumx += xx * dum;
X            }
X        }
X
X        stemp = (sumx / sum) - a;
X        if ((np-nm) >= 2) {
X            /* guess is too low make another pass */
X            am = a;
X            aa = stemp < 0.0? xp : xp + stemp*AMP;
X            if (aa > ap) aa = 0.5*(a+ap);
X            eps = AFAC * fabs((double)(aa-a));
X            a = aa;
X        } else if ((nm-np) >= 2) {
X            /* guess is too high make another pass */
X            ap = a;
X            aa = stemp > 0.0? xm : xm+stemp*AMP;
X            if (aa < am) aa = 0.5*(a+am);
X            eps = AFAC*fabs((double)(aa-a));
X            a = aa;
X        } else {    /* got it    */
X            if (NTuples%2 == 0) {
X                Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
X            } else {
X                Median = np == nm? a : np>nm? xp : xm;
X            }
X            return(Median);
X        }
X    }
X}
X
XFLOAT
XMedianOfY(Data,N,M)
XFLOAT    **Data;
Xint    N;
Xint    M;
X/* an iterative computation of the median...
X * Code taken from Numerical Recipes..
X */
X{
X    int    n, np, nm, i,j,D;
X    FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
X    FLOAT    Median;
X    double    fabs();
X
X    D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
X    n = N * (M-D);
X
X    /* first estimate    */
X    a = 0.5 * (Data[0][D] + Data[N-1][M-1]);
X    eps = fabs((double)(Data[N-1][M-1] - Data[0][D]));
X
X    am = -(ap=BIG);
X    for (;;) {
X        sum = sumx = 0.0;
X        np = nm = 0;        /* # point above and below median */
X        xm = -(xp = BIG);
X
X        for (i=0; i<N; i++) {
X            for (j=D; j<M; j++) {
X                xx = Data[i][j];
X                if (xx != a) {
X                    if (xx > a) {
X                        ++np;
X                        if (xx < xp) xp = xx;
X                    } else if (xx < a) {
X                        ++nm;
X                        if (xx > xm) xm = xx;
X                    }
X                    sum += dum=1.0/(eps+fabs((double)(xx-a)));
X                    sumx += xx * dum;
X                }
X            }
X        }
X
X        stemp = (sumx / sum) - a;
X        if ((np-nm) >= 2) {
X            /* guess is too low make another pass */
X            am = a;
X            aa = stemp < 0.0? xp : xp + stemp*AMP;
X            if (aa > ap) aa = 0.5*(a+ap);
X            eps = AFAC * fabs((double)(aa-a));
X            a = aa;
X        } else if ((nm-np) >= 2) {
X            /* guess is too high make another pass */
X            ap = a;
X            aa = stemp > 0.0? xm : xm+stemp*AMP;
X            if (aa < am) aa = 0.5*(a+am);
X            eps = AFAC*fabs((double)(aa-a));
X            a = aa;
X        } else {    /* got it    */
X            if (n%2 == 0) {
X                Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
X            } else {
X                Median = np == nm? a : np>nm? xp : xm;
X            }
X            return(Median);
X        }
X    }
X}
END_OF_FILE
if test 25452 -ne `wc -c <'Csrc/plotting.c'`; then
    echo shar: \"'Csrc/plotting.c'\" unpacked with wrong size!
fi
chmod +x 'Csrc/plotting.c'
# end of 'Csrc/plotting.c'
fi
echo shar: End of archive 3 \(of 3\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 3 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0
--
Mail submissions (sources or binaries) to <[e-mail]>.
Mail comments to the moderator at <[e-mail]>.
Post requests for sources, and general discussion to comp.sys.amiga.